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Abstract 

We present a theory of electron, electromagnetic, and elastic wave propagation in 
systems consisting of non-overlapping scatterers in a host medium. The theory pro- 
vides a framework for a unified description of wave propagation in three-dimensional 
periodic structures, finite slabs of layered structures, and systems with impurities: 
isolated impurities, impurity aggregates, or randomly distributed impurities. We 
point out the similarities and differences between the different cases considered, 
and discuss the numerical implementation of the formalism. 
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1 Introduction 



The Korringa-Kohn-Rostoker (KKR) method [1] has been used extensively in 
the study of the electronic structure and related properties of materials. It 
has been very successful in calculations of the electronic structure of ordered 
elemental solids [2], impurities [3], and disordered alloys [4]. Based on the 
multiple-scattering theory [5], the KKR method calculates from first princi- 
ples quantities such as the electron density and the ground-state energy of the 
system. A KKR method for layered systems, the so-called layer KKR (LKKR) 
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method, has also been developed [6,7]. This method is well adapted to calcu- 
lations of low-energy electron diffraction (LEED) [8] and electron emission [9] 
spectra. 

In recent years, propagation of classical waves [electromagnetic (EM) or elastic 
waves] in composite materials with dielectric or, respectively, elastic proper- 
ties which are periodic functions of the position, with a period comparable to 
the wavelength of the corresponding field, has been the object of considerable 
attention [10]. These materials, photonic and phononic crystals, respectively, 
whether they exist naturally or are artificially fabricated, exhibit a rich variety 
of physical properties of interest to fundamental and applied research. There 
are striking analogies between the propagation of electrons in ordinary crystals 
and EM/elastic waves in photonic/phononic crystals, so that a great variety 
of multiple-scattering methods originally developed for electronic-structure 
calculations have been transferred to the field of photonic and phononic crys- 
tals [11-15]. The aim of the present article is to present the KKR formalism 
for electron, EM and elastic waves from a unifying point of view. 



2 Multipole expansion of a wave field 



We first consider, as an example of a field with one degree of freedom, the case 
of a spinless particle of mass m and energy £, in a potential V{r). According 
to quantum mechanics, the particle is described by a scalar field: the single- 
particle wave function 

^(r; t) = ^(r) cxp {-iSt/h) , (1) 



where k is the Planck constant and ^(r) satisfies the time-independent Schrodinger 
equation [16] 



* = 0. 



(2) 



In a region of constant potential, Vq, Eq. (2) reduces to the Helmholtz wave 
equation 



+ g^l ^ = , 



(3) 



where a = ^j2m{E -Vo)/h. A complete set of spherical-wave solutions of 
Eq. (3) is given by 



^,m{r) = fe{qr)Yr{r), 



(4) 
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where Yf^{r) are the usual spherical harmonics and fc may be any linear 
combination of the spherical Bessel function, jf, and the spherical Hankel 
function, h^. The most general wave function in a constant potential field can 
be written as hnear sum of the spherical waves given by Eq. (4), as follows: 

*(r) = $:a,„/,(gr)rr(f), (5) 

im 

where a^m are coefficients to be determined. 

A vector field has in general, not one, but three degrees of freedom. Let us, 
for example, consider a harmonic elastic wave, of angular frequency lu, propa- 
gating in a homogeneous medium of mass density p and Lame coefficients A, 
fjL. This is described by a displacement vector field: 

U(r, t) = Re [U(r)exp(-ic;t)] , (6) 

where U(r) satisfies the equation [17] 

(A + 2//) V (V • U) - AtV X (V X U) + fxj'^U = . (7) 



The most general solution of Eq. (7) consists of two types of elastic waves 
which propagate independently. These arc purely longitudinal (irrotational) 
waves (V X U = 0), which satisfy the vector Helmholtz equation 



V' + qf 



U = 0, 



(8) 



where qi = lo/q, with q = y(A + 2yu)/p being the longitudinal-wave speed; 
and purely transverse (divergence-free) waves (V • U = 0), which satisfy the 
vector Helmholtz equation 



(9) 



where qt — cu/ct, with q = y/jjjp being the transverse- wave speed. A complete 
set of longitudinal (L) spherical- wave solutions of Eq. (8) is given by [18] 



UL(r) 



Qi 



-V[/,(g;r)r/ 



(10) 



A complete set of transverse (M, N) spherical-wave solutions of Eq. (9) is 
given by [18] 



U,^(r) = fe{qtr)Xem{r) 
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and 

UL(r) = -V X /,(g,r)X,^(f) , (12) 



where X.£rn{^) are the so-called vector spherical harmonics, which are defined 

by 

^£(£ + l)X,^(f) = Lrr(f) = -ir x Vrr(r) ■ (13) 

The most general displacement field in a homogeneous medium can be written 
as a linear sum of the spherical waves given by Eqs. (10)-(12), as follows 

U(r) = J2{^ZrMQtr)^im{f) + af„-V x /,(gir)X,^(f) 
em ^ 

+aL-V[/,(ftr)Fr(f)]j, (14) 

where a^, P — M,N, L, are coefficients to be determined. 

Another example of a vector field is the EM field. A harmonic EM wave, of 
angular frequency a;, is described by its electric-field component 

E(r, t) = Re [E(r)exp(-ia;t)] . (15) 

We need not write down explicitly the magnetic-field component of the wave 
since this can be readily obtained from E(r, t) using one of Maxwell equations. 
In a homogeneous medium characterized by a dielectric function e(a;)eo and a 
magnetic permeability fi{uj)fio, where Cq, /io are the electric permittivity and 
magnetic permeability of vacuum, Maxwell equations imply that E(r) satisfies 
a vector Helmholtz equation, subject to the condition V ■ E = 0, with a wave 
number q — uj/c, where c = 1 / ^//Ie/7oeo is the velocity of light in the medium. 
The spherical- wave expansion of E(r) is given by the terms of Eq. (14) which 
satisfy V • E = 0, i.e., 

E(r) = E |af„/,(gr)X,^(r) + a.^-V x [/,(gr)X,^(f)]| , (16) 
em I ^ ) 

where a^, P = H,E, are coefficients to be determined. 
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3 Scattering by a single scatterer 

We consider a scatterer of finite range, S, with its center at the origin of coordi- 
nates, and assume that the appropriate in each case characteristics (potential. 
Lame coefficients and mass density, electric permittivity and magnetic perme- 
ability) are different from those of the surrounding homogeneous medium. A 
De Broglie, EM or an clastic plane wave incident on this scatterer is described, 
respectively, by Eq. (4), (16) or (14) with fe — (since the plane wave is finite 
everywhere) and appropriate coefficients a° , where L denotes collectively the 
indices imP. Similarly, the scattered wave is described by Eq. (4), (16) or 
(14) with fe = h'l, which has the asymptotic form appropriate to an outgoing 
spherical wave: hj^ ~ (— i)^exp(igr)/igr as r — > cxd, and appropriate expansion 
coefficients aj. The wavefield for r > is the sum of the incident and scat- 
tered waves. The spherical-wave expansion of the field for r < yS" is obtained 
in similar manner by the requirement that it be finite at the origin (r = 0). 
By applying the proper in each case boundary conditions at the surface of 
the scatterer, we obtain a relation between the expansion coefficients of the 
incident and the scattered field, as follows: 

L' 

where Tll' are the elements of the so-called scattering transition matrix. 
Eq. (17) is valid for any shape of scatterer; for spherically symmetric scat- 
terers each spherical wave scatters independently of all others, which leads to 
a transition matrix diagonal in angular momentum. Explicit forms for Tl^i 
for various cases can be found elsewhere [16,19,20]. 



4 Multiple scattering theory 

We consider an assembly of non-overlapping scatterers centered on sites Rj 
in a homogeneous medium. In general: an outgoing wave about [described 
by Eq. (4), (14) or (16) with fi = hj' and expansion coefficients b^* ] can be 
expanded into spherical waves about Rj, incident on Rj. This expansion has 
the form of Eq. (4), (14) or (16) with = ji and expansion coefficients b'^(i'), 
given by 

b'l(i')=E^li'bir. (18) 

L' 

Explicit expressions for the so-called free-space structural Green functions, 
for the scalar and transverse vector fields can be found elsewhere [4,21]. 
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In the case of the elastic field, since a longitudinal (transverse) spherical wave 
about Rj/ remains a longitudinal (transverse) wave when expanded about 
another center Rj 7^ Rj/, the matrix elements ^lli are obtained independently 
for longitudinal and transverse waves. Therefore, Q^l' ^ block-diagonal 
form in the mode index P (we recall that L = £mP): the matrix elements of 
the block corresponding to transverse waves are identical to the corresponding 
matrix elements for the EM field, while the matrix elements of the block 
corresponding to longitudinal waves are identical to the corresponding matrix 
elements for the scalar field [15]. We note that, for any field, by definition, 
^LL' equals zero for i — i'. 

The wave scattered from the scatterer at Rj is determined from the total wave 
incident on this scatterer; therefore 

hr^T.TLL'i^'^' + <'), (19) 

L' 



where a°' are the coefficients in the multipole expansion [given by Eq. (4), (14) 
or (16)] about R^ of an externally incident wave; and b'^ are the coefficients 
in the multipole expansion about Rj of the outgoing waves scattered from all 
the other scatterers at Rj/ 7^ Rj. From Eq. (18) we have 



b'V 



E 

i',L' 



Qii' Y^+i 



(20) 



Substituting Eq. (20) into Eq. (19) we obtain 



E U'^LL' - E tIl" K^' = E Ul' ■ (21) 

i',L' V L" J L' 



We now introduce the structural Green functions -D|,£,/ which give the coef- 
ficients [in an expansion such as (4), (14) or (16)] of the wave incident on 
the scatterer at Rj, due to an outgoing wave from the scatterer at Rj/. An 
outgoing wave from the i'th scatterer can reach the ith scatterer directly, or 
indirectly after scattering any number of times by any number of scatterers 
(including those at Rj and Rj/). Let express the sum of the contribution 
to the coefficients of the incident wave on Rj from all possible scattering paths 
originating from an outgoing wave from the i'th scatterer. One can easily prove 
by iteration the following equation 

^LV — ^LL' + E ^LL"'^L" L'" -^L'" L' i (^2) 
i",L",L"' 



6 



where Tl'',j^,„ are the elements of the transition matrix which describes the 
scattering by the single scatterer at Rj". We generalize Eq. (22) by treating 
the scattering at Rj in two stages. The first-stage scattering is described by 
T£^,, which correspond to arbitrarily defined scatterers (reference scatterers), 
and the second-stage scattering by T/^, — T£^,. We obtain: 

T-jii' 7-)ru' I \ ^ T^rii" (rni" rpri" \ T^i"i' (OQ\ 

^LL' — ^LL' "T ^LL" y'- L"L"' L"L"' ) ^L"'L' ) K'^'^) 

i",L",L'" 



where i^^^/ are the solution of Eq. (22) when T£x,/ = Tj 



5 Bulk systems 



We consider the case in which the scatterers constitute a three dimensional 
(3D) crystal structure specified by Bravais lattice vectors Rn and non-primitive 
translation vectors, indicating the positions of the scatterers within the unit 
cell, t„; in this case the site index i stands for the composite index na. Bloch's 



theorem implies that b^^ " = exp [ik 



(R„/ — R^ 



The normal modes of the crystal are obtained, by putting the external incident 
wave equal to zero in Eq. (21), which leads to the following secular equation: 



det 



L" 



0, 



(24) 



where 

^II' (k) = E exp [ik • (R„, - R„)] , (25) 

n' 

which does not depend on n. Both T£^, and (k) in Eq. (24) are func- 
tions of the frequency of the wave, but the T^"^, depend only on the properties 
of a single scatterer, whereas Q"£/(k) depend only on the geometry. Thus 
Eq. (24) reflects a complete separation of the individual scatterer and struc- 
tural aspects of the problem. The matrix elements (k) for the scalar field 
were introduced by Korringa, Kohn, and Rostoker [1] in relation to the cal- 
culation of the electronic band structure of periodic solids and are commonly 
referred to as the KKR structure constants. The calculation of the structure 
constants, which needs to be done only once for a given lattice, usually re- 
quires Ewald-summation techniques [22]. We note that Eq. (24) involves in 
principle infinite-dimensional matrices. In actual calculations, however, it is 
sufficient to truncate the angular momentum index, £, to some relatively small 
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number, imax- The KKR method has been successfully used for many years in 
the calculation of the electron band structure of solids [2], and more recently 
it has also been applied to the calculation of the frequency band structure of 
EM and elastic fields in relation to photonic and phononic crystals [14]. 

One of the advantages of the KKR theory compared to other approaches is 
its capability of dealing with defects and disorder [4]. For the description of 
point defects at a finite number of sites, one needs to calculate the structural 
Green functions, Dj^j^,, of the defect system. This can be done in real space, 
using Eq. (23) and considering the periodic crystal (without any defect) as 
the reference system. In this case, the sum over i" in Eq. (23) is restricted 
to those sites at which there are defects and induce a perturbation STl„j^,„ = 
Tl'iL'" ~ Tyiyi, on the T matrix. For the periodic arrangement of scatterers 
of the reference system, the evaluation of D^l' through Eq. (22) involves an 
infinite number of sites i" and it is achieved by a lattice Fourier transform of 
Eq. (22) 

Dll?\\,) = nti{\,)+ f^rL'"(k)TI,^'..i^I^:i^'(k) (26) 

a",L",L"' 



and subsequent integration over the volume v of the Brillouin zone (BZ) 

D^, = - [d'k exp [ik ■ (R„ - R„0] Dlt?'{k) . (27) 
V J 
BZ 



If a number of N different scatterers, described by T matrices Tl^,, j ~ 
1, 2, . . . , A^, are randomly distributed on the sites of a given crystal, the coher- 
ent potential approximation (CPA) [23] can be used to treat the disordered 
system within a mean- field context. The effective CPA medium consists of 
identical scatterers at all sites of sublattice a of the crystal, characterized 
by a T matrix Tlf^, . The structural Green functions D'l^"'"' " of the CPA 
medium are calculated by equations analogous to (26) and (27). The CPA 
self-consistency condition demands that the correction to the scattering due 
to the difference of the actual scatterers from the CPA scatterer vanishes on 
the average 



^ 1-1 



Cj-a [l - (t^' - T'^") °";°"] (t^' - T'^") = for every a , (28) 



where Cj.^a denotes the concentration of scatterer j on sublattice a. 
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6 Layered systems 



A layered system consists of a number of planes (layers) of scatterers with the 
same two dimensional (2D) periodicity. To begin with, we consider just one 
layer, at 2; = 0, in which case the scatterers are centered on the sites R„ + 1^, 
where {Rn} constitutes a 2D Bravais lattice and t„ denote the positions of 
the scatterers within the 2D unit cell. The 2D reciprocal vectors g, and the 
surface Brillouin zone (SBZ) corresponding to this lattice are defined in the 
usual manner [8,9]. The wave vector parallel to the plane of atoms can always 
be written as follows: 

q||=k||+g', (29) 



where the reduced wave vector, ky , hes in the SBZ and g' is a certain reciprocal 
vector. In what follows we write the wave vector of a plane wave of given 
qy = k|| + g and given wave number q^,, where v specifies the polarization 
mode of the wave, if there exist longitudinal and transverse modes, as follows: 



K 



± 



1/2 



(30) 



We note that when < (ky + g)^, the above defines a decaying wave; the 
positive (negative) sign in Eq. (30) corresponds to a wave propagating or 
decaying to the right (left). 



The coefficients b^™"^ in the multipole expansion of the scattered field from 
the layer are obtained from Eq. (21). Using the Bloch condition: b^^" = 
exp ik|| • (Rjj/ — R„) b^"" , we obtain 



E 

a',L' 



L 



,Qna 



(31) 



where a^^"" 



wave. 



are the coefficients in the multipole expansion of the incident plane 
The 2D KKR structure constants 



(k||)=E^ 



na\n a 
LL> 



exp 



ikii • (R„/ — Rt, 



(32) 



hke their counterparts in 3D given by Eq. (25), can be evaluated using Ewald- 
summation techniques [6]. 

Writing the incident, refiected and transmitted plane waves with respect to 
an origin (a scattering center) in this plane, we obtain the amplitudes, [$rf]gi 
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and [$tr]gj of the reflected and transmitted beams, respectively, in terms of 
the amphtudes, [$in]g/j, of the incident wave as follows: 

i' 

[*tr]^=E^gl;gv[*in]g'.- (33) 

i' 

where s = +(— ) specifies a wave travehng or decaying to the positive (neg- 
ative) z direction, and i denotes the components of the wave field. A scalar 
field — > ^) has just one component {i — 1), the displacement vector 
of the elastic field — > U) has three components {i — 1,2,3), while the 
EM field — > E), because of its transverse nature, has only two indepen- 
dent components, i = 1,2. Exphcit expressions for the transmission/reflection 
matrix elements M||.g,j, for the scalar, elastic, and EM fields can be found 
elsewhere [8,9,15,24].' 

The transmission/reflection matrices for a slab which consists of a stack of 
layers with the same 2D periodicity parallel to a given plane are obtained from 
the transmission/refiection matrices of the individual layers in the manner 
described in Refs. [8,9,15,24]. Knowing the transmission/refiection matrices 
for the slab we can readily obtain the transmission, reflection, and absorption 
coefficients of the plane wave of Eq. (30) incident on the slab. The LKKR 
technique summarized in this section is particularly well adapted to layered 
systems because the computation time scales linearly to the number of layers 
in the slab. Moreover, because this method does not require periodicity in the 
direction perpendicular to the layers, it can easily treat heterostructures, slabs 
with impurity planes [24-27] or even a random succession of different layers, 
as long as these have the same 2D periodicity. 

The LKKR technique provides also the complex band structure of an inflnite 
crystal, viewed as a succession of layers parallel to a given crystallographic 
plane [8,9,15,24]. 

Point defects or a random distribution of point defects in the layers of a slab 
can also be treated along the lines described in Section 5 [4,21,28]. 
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